Elasticity of Filamentous Kagome Lattice 
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The diluted kagome lattice, in which bonds are randomly removed with probability 1 — p, consists 
of straight lines that intersect at points with a maximum coordination number of four. If lines are 
treated as semi-fiexible polymers and crossing points are treated as crosslinks, this lattice provides 
a simple model for two-dimensional filamentous networks. Lattice-based effective medium theories 
and numerical simulations for filaments modeled as elastic rods, with stretching modulus /i and 
bending modulus k, are used to study the elasticity of this lattice as functions of p and re. At p — 1, 
elastic response is purely afline, and the macroscopic elastic modulus G is independent of re. When 
k = 0, the lattice undergoes a first-order rigidity percolation transition at p — 1. When re > 0, G 
decreases continuously as p decreases below one, reaching zero at a continuous rigidity percolation 
transition at p — pb ~ 0.605 that is the same for all non-zero values of re. The effective medium 
theories predict scaling forms for G, which exhibit crossover from bending dominated response at 
small re//i to stretching-dominated response at large re//i near both p = 1 and p — pb, that match 
simulations with no adjustable parameters near p = 1. The afline response as p — > 1 is identified 
with the approach to a state with sample-crossing straight filaments treated as elastic rods. 

PACS numbers: 87.16. Ka, 61.43.-j, 62.20.de, 05.70.Jk 



I. INTRODUCTION 

Filamentous networks [1] are an important class of 
materials in which the interplay between stretching and 
bending energies and temperature lead to unique me- 
chanical properties such as strong nonlinear response and 
strain stiffening [2-5], negative normal stress [5, 6], non- 
affine response [4, 7-9], crossover from non-affine to affine 
response [10, 11], and power-law frequency dependence 
of the storage and loss moduli [12]. They are a part 
of such important components of living matter [13-16] 
as the cytoskeleton, the intercellular matrix, and clotted 
blood and of industrial materials like paper [17, 18]. Here 
we introduce the diluted kagome lattice, shown in Fig. 1, 
in which elastic rods on nearest-neighbor (NN) bonds 
are removed with probability 1 — p, as a model for fila- 
mentous networks in two-dimensions. Lines of contiguous 
and colinear occupied bonds are identified as filaments, 
which arc modeled, as in previous work [7-11, 19] as elas- 
tic rods with one-dimensional stretching modulus /i and 
bending modulus k. 

Effective medium theories (EMTs) [20-22] have proven 
to be a powerful tool for the study of random systems. 
We use recently derived EMTs [23-26] that treat bending 
as well as stretching forces in elastic networks to calcu- 
late the shear modulus G as a function of /it, re, and p, 
and derive scaling forms for its behavior near p = 1 and 
near the re > rigidity threshold at p — pi,. We also cal- 
culate G using numerical simulations. The simulations 
and EMTs are in qualitative agreement over the entire 
range of p and in quantitative agreement near p = 1. 
Our results are also in general agreement with those for 
the Mikado model [4, 10, 19], including, in particular, a 
crossover from bending-dominated nonaffine response at 
low density (small p) to stretching dominated affine re- 




FIG. 1. A bond-diluted kagome lattice with p — 0.8 and 
N x = 14 horizontal bonds and N y = 12 bonds along the axis 
at an angle a = n/3 to the x-axis. The red sites on the lattice 
are crosslinks, and colinear sequences of bonds are filaments. 
Only two filaments, shown with thicker bonds, traverse the 
lattice from top to bottom and support external stresses. One 
these makes an angle 4>\ = n /3 = a and one an angle <j>2 = 
2-7r/3 = 2a with the a;-axis, and the length of both is N y a. In 
more general lattices, sample crossing filaments can orient at 
any angle <j> to the x-axis, but the projection of their length 
onto the vertical axis is W y , and thus their length is L(<j>) = 
W y /{a sin^). 



sponse at high density (p near 1). Our model, however, 
highlights special properties of networks of straight fila- 
ments, including scaling collapse, described analytically 
by our EMTs, of the region near p = 1 where filament 
length L approaches infinity, and the underlying origins 
of affine response in this limit. 

The building blocks of filamentous networks are typi- 
cally semi-flexible polymers, characterized by a bending 
modulus k = IpksT, where l p is the persistence length kg 
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the Boltzmann constant, and T the temperature. Where 
filaments intersect, physical or chemical crosslinks bind 
them together but do not inhibit their relative rotation. 
Since only two filaments cross at a point, each crosslink 
is connected to at most four others. The average distance 
l c between crosslinks is less that l p . A remarkably faith- 
ful description of experimentally measured linear [27, 28] 
and nonlinear [3] elastic response is provided by a model 
in which elastic response is assumed to be affine and in 
which filamentous sections between crosslinks are treated 
as nonlinear central-force springs with a force-extension 
curve determined by the worm-like chain model [3, 27, 29] 
augmented by an enthalpic stretching energy. This sim- 
ple model, however, misses some important properties of 
filamentous networks, including the existence of a rigidity 
percolation threshold [30, 31], nonaffine response [4, 7- 
9], and crossover from bending-dominated nonaffine re- 
sponse at small L/l c [10, 11], or low frequency [32], to al- 
most affine, stretching dominated response at large L/l c 
or high frequency. 

The Mikado model [11, 19] provides additional insights 
into the physics of filamentous networks. In this model, 
straight filaments of length L are deposited with random 
position and orientation on a two-dimensional flat sur- 
face and crosslinkcd at their points of intersection. Fil- 
aments are treated as elastic rods with one-dimensional 
stretching and bending moduli, /x and k, rather than as 
worm-like chains. Numerical simulations on this model 
reveal a rigidity percolation transition from a floppy net- 
work to one with non-vanishing shear and bulk moduli 
and with non-affine bending dominated elastic response. 
As L/l c is increased, response becomes more affine and 
stretching dominated, reaching almost perfectly affine re- 
sponse in the L/l c — > oo limit. Our EMTs and numerical 
simulations yield similar results for the diluted kagome 
model. 

Though the Mikado and the diluted kagome lattice 
are quite similar with crosslinks of maximum coordi- 
nation four and variable values of the ratio L/l c , they 
model slightly different parts of the phase space of pos- 
sible two-dimensional filamentous networks. In particu- 
lar, since the starting point of the diluted kagome lat- 
tice is the full lattice with all bonds occupied, it nec- 
essarily deals directly with the L — > oo limit, which is 
not generally accessed in studies of the Mikado model in 
which L is restricted to be less than the sample width W . 
The Mikado model is an "off-lattice" model, whereas the 
kagome model is lattice based. The latter property of the 
kagome model facilitates the application of lattice-based 
EMTs [20-22] (modified to include bending) for all values 
of p. Finally, In the Mikado model, the distance between 
between crosslinks on a single filament follows a Poisson 
distribution with no lower bound, whereas in the kagome 
lattice, this distance is an integral multiple of the lattice 
spacing with a lower bound of one lattice spacing. This 
difference leads, as we shall see, to a different scaling of 
the shear modulus in the bending dominated regime near 
p = l. 



As detailed in Apps. B and C, our EMTs clearly show 
that there are three critical points (or fixed points in the 
rcnormalization-group sense): the trivial empty lattice 
point at p — (which we ignore) , the rigidity-percolation 
point at p = pi,, and the p = 1 full lattice point. It 
also provides analytic scaling relations for the the shear 
modulus G as a function of fi, k, and p in the W — > oo 
limit near both p — Pb and p = 1. Of particular interest 
is the behavior near p = 1, where the ratio G/Go of the 
shear modulus to its p = 1 value Go can be expressed as 
a function of the single scaling variable 



T = 



(1 — p) 2 fia 2 



(i.i) 



where l\ = njp, is the bending length (L) = a/(l — p) 
is the average filament length and l c w a is the average 
spacing between crosslinks along a filament. As r — > oo, 
G approaches Go, but as r — > 0, there is crossover to a 
bending dominated, nonaffine regime in which 



G 



(1.2) 



This behavior, which appears in our simulations, has also 
been observed in simulations in diluted three-dimensional 
lattices of four-fold coordinated filaments [33, 34]. In 
Sec. V, we speculate about the reasons for the same 
behavior appearing in both two- and three-dimensions. 
Huessinger et al. [7, 8] developed an off-lattice EMT for 
filamentous networks that predicts a bending-dominated 
non-affine regime with a different power of L/l c , G ~ 
(k/1^,){L/I c ) a , than the one we predict. The origin of 
this different scaling is the absence of a lower cutoff in the 
Poisson distribution of the distance between crosslinks in 
the Mikado model compared to the fixed cutoff equal to 
the lattice spacing in the kagome model. The analysis 
of Refs. [7, 8] yields the kagome lattice scaling law if the 
probability distribution for distances is replaced by one 
with a fixed lower cutoff. 

Our numerical simulations follow the EMT prediction 
very closely with no adjustable parameters near p = 1 
if k is not too large. For really small values of k//j,o, 2 , 
finite-size effects become important, and simulations can 
be fit with a combination of the exactly calculable finite 
size results at k — and the EMT scaling form valid 
at infinite W. Near p — pb, simulations are consistent 
with a scaling function of the form of Eq. (??) but with t 
approximately 0.2 rather than 1 and with different values 
of ci and c 2 . 

Under periodic boundary conditions, the kagome lat- 
tice has exactly z = 4 = 2d neighbors per site, where d is 
the spatial dimension. If sites interact only via central- 
force springs on bonds and not by bending forces along 
filaments, the lattice is isostatic in that it has exactly 
enough internal forces to eliminate zero modes accord- 
ing the simple Maxwell count [35]: Nq — dN — Nb = 
(d — 7}Z)N, where Nq is the number of zero modes and 
Nb = \zN is the number of bonds. Under free rather 
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FIG. 2. Schematic phase diagram for the diluted kagome lat- 
tice showing the bending-dominated region in which G ~ « 
and the stretching dominated region in which G ~ /i. This 
boundary between stretching- and bending-dominated re- 
sponse was calculated using a formula that interpolates be- 
tween the analytic EMT expressions for G near p = pb and 
near p = 1. 



than periodic boundary conditions, sites on the bound- 
aries have fewer than four neighbors, and there is an over- 
all bond deficiency of AN B = IN - N B ~ VN. As a re- 
sult, there are of order VA zero "floppy" internal modes 
of zero energy in the absence of bending forces. Under 
periodic boundary conditions, the simple Maxwell count 
leads to A = 0, but the total number of zero modes 
(as calculated by the dynamical matrix, whose normal 
modes are the "infinitesimal" modes that result when 
the elastic energy is truncated to harmonic order), like 
the number under free boundary conditions, is of order 
y/~N rather than zero [36]. This discrepancy is a result 
of the existence of lattice configurations that can sup- 
port stress in bonds while maintaining a zero force at 
each note [37]. Associated with each such configuration 
[38], called a "state of self stress" , is an additional in- 
finitesimal zero mode so that the Maxwell count becomes 
N n = oIN — Nb + S, where S is the number of states of self 
stress. Thus, in the undiluted central-force kagome lat- 
tice under periodic boundary conditions, Nq — S is sim- 
ply the number of states of self-stress. These states en- 
dow the undiluted central-force kagome lattice with affine 
response and non-zero bulk and shear moduli, propor- 
tional to the bond spring constant, and, thus, determine 
the approach to affine, stretching-dominated response as 
p — > 1 in the diluted lattice with nonzero k. In addition, 
the central-force zero modes associated with the states of 
self stress also control the form of the effective-medium 
equations and the scaling forms they predict. The ad- 
dition of bending forces lifts all but the two trivial zero 
modes of rigid translation to finite frequency. This fact 
plays a central roll in our EMT and is ultimately respon- 
sible for the form of the scaling function for G near p = 1. 

We are ultimately interested in three-dimensional fil- 
amentous networks, which are subisostatic with Nq = 
3 A — 2N = N because their maximum coordination num- 
ber, like the that of two-dimensional networks, is four. If 



constituent filaments are straight, these networks have a 
state of self-stress for each filament in the undiluted limit, 
and their elastic moduli are non-zero, proportional to the 
bond-spring constant, and independent of K [33]. Associ- 
ated with each state of self-stress, there is an additional 
zero mode, but as in the kagome lattice, the addition of 
bending forces raises all but the three trivial zero modes 
of uniform translation to finite frequency and stabilizes 
the lattice. We will argue in Sec. V that these properties 
are the likely origin of the very similar behavior, seen in 
simulations, of the shear modulus in 2- and 3-dimensional 
networks with maximum coordination of four. 

This paper consists of five sections and three appen- 
dices. Section II compares and contrasts the Mikado and 
kagome models and demonstrates how straight, sample- 
traversing filaments with the energy of an elastic beam 
produce affine elastic response. Section II also defines 
the lattice model that we use. Section III describes 
EMT procedures and presents their results. Section IV 
presents the results of our numerical simulations and 
compares them with those of the EMTs. Section V re- 
views our results and speculates about application of our 
two-dimensional calculations to three-dimensional sys- 
tems. There are three appendices: Appendix A discusses 
general properties of the EMT dynamical matrix. Ap- 
pendices B and C present details of the solutions to the 
EMT equations near p — 1 and p = pt for the bending 
EMTs described in Refs. [24, 25] and in Refs. [23, 26], 
respectively. 

II. FILAMENTOUS NETWORKS AND THE 
DILUTED KAGOME LATTICE 

Networks composed of straight filaments with elastic- 
rod energies have special properties. In this section, we 
will explore some general properties of these networks 
before setting up the energy for the kagome lattice. 

A. Elastic filamentous networks 

The elastic energy of a one-dimensional elastic beam 
with stretching modulus fj, and bending modulus k is 

(2.1) 

where u(s) is the local longitudinal displacement and 9(s) 
the local angle of the tangent to the filament at the point 
s. The usual practice is to treat the stretching modulus as 
independent of k even though at nonzero temperature the 
spring constant for the entire filament has an important 
entropic component that depends on k . 

The spring constant for a filament section of length I is 
k = fil^ 1 . This leads to affine response in linear elasticity 
for any lattice in any dimension consisting of sample- 
traversing straight filaments with a sufficient number of 
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orientations to ensure stability with respect to all strain 
[7, 8, 39]. To see this, consider a crosslink (node of the 
lattice) at the origin on a filament parallel to the unit 
vector e, and let it be connected to two other crosslinks 
on the same filament at respective positions Si = lie and 
si = — Z 2 e relative to the orign. Under a uniform, affine 
deformation, the relative positions s a , a — 1,2 transform 
to A • s a = s a + rj ■ s a , where A = / + 77 is the deformation 
tensor with components Ay = 5ij + 77^- . The forces that 
crosslinks 1 and 2 exert on the origin are then 

fi = k(h)(heir]ijej)e; f 2 = -k(fo)(foeirnjej)e, (2.2) 

where i and j run over x, y and the summation conven- 
tion is understood. The sum of these forces is zero be- 
cause k{l) ~ n/l. The same analysis applies to any site 
and filament in the lattice. Under affine distortions, fila- 
ments do not bend, so the energy of the affine distortion 
depends only on the central force and does not depend 
on k. Thus, under affine distortions of sample travers- 
ing filaments, the force on every intermediate crosslink is 
zero, and non-affine distortions are not introduced: the 
response is affine and independent k. When the lattice 
is diluted, not every crosslink is connected to two others, 
the above cancelation does not occur, and the result is 
nonaffine response with a bending component. 

With these observations, we can calculate the shear 
modulus of any network of straight filaments with 
stretching energy only, i.e., k = 0. For simplicity, we 
restrict our attention to two-dimensional systems in a 
rectangular- or rhombus-shaped box with base-length W x 
and height W y as shown in Fig. 1, and we consider only 
shear deformations in which the only nonvanishing com- 
ponent of 77 is r] xy . In this case, only those filaments 
that cross from the bottom to the top of the sample will 
contribute to the shear modulus. The length of such a 
sample traversing filament oriented along a unit vector 
e(0) = (cos <f>, sin <fi) making an angle <p with the x-axis 
is L(<j>) = W y I '(a sin cf>). Under a shear deformation 77, its 
length will change by 5L = Lei((j>)riijej((f)) to lowest or- 
der in 77. Since the spring constant of a filament of length 
L is /i/L, the clastic energy of the filament is 

EM = 1 (2.3) 

_ bill G) 

and the total energy from all filaments is E to t — 
(1 / '2) NW y ii{\ei{4>)r]ij ej (<f))} 2 / sin <f>), where N is the total 
number of sample traversing filaments and (• • • ) signifies 
an average over the orientation angles of the filaments. 
The linearized energy density, E/W x W y , is 

/ = \ Ci jk i UijU U , (2.4) 

where Cijki are the components of the elastic constant or 
elastic modulus tensor and Uij — are the com- 

ponents of the linearized strain tensor u. In the isotropic 
limit, the energy density reduces to 

/ = ^B(Tr^) 2 + GTr^ 2 , (2.5) 



where u is the traceless part of u, and G = C xyxy and 
B = \{C XXXX + C xxyy ) are the shear and bulk moduli, 
respectively. These relations along with Eq. (2.3) for the 
case in which the only nonvanishing component of rjij is 
rj xy lead to 

G = I^v/!^\ (2 . 6) 

AW X \ sin0 / k ' 

for the shear modulus of a network of elastic rods with 
stretching modulus /j. Note that rods parallel to the x- 
axis {(j) — 0) do not contribute to G. In the kagome lat- 
tice, bottom-to-top traversing filaments all have length 
L v = W y / sin(7r/3) and their angles with respect to 
the x-axis are restricted to <fi = 7r/3, 27t/3 for which 
sin</> = sin 20 = V3/2. There are W x /a sites along the x 
axis from which a single filament aligned cither along n/3 
or 27r/3 can emerge, so N is simply g cross W x /a, where 
Across ip) is the probability a filament emerging from a 
given site along the x-axis traverses the sample. For sim- 
plicity, we take L y = W in which g cr0 ss(p, W) = p w / a . 
When p = 1, all bonds are occupied, and the shear modu- 
lus of the undiluted kagome lattice is Go = (\/3/8)(/u/a). 
When p < 1 and k = 0, there is a nonzero stretching 
contribution to G in finite samples: 

G s tA^P, W) = G(h,k = 0,p, W) = q CIOSS (p, W)G . 

(2-7) 

In the W — > 00 limit, the probability of sample travers- 
ing filaments vanishes for any p < 1, and Goo(/J, n,p) — 
limv^-j-oo G(fi, K,p, W) is zero at k = for all p < 1. 
Thus, — > 0, G rx (fi, K,p) must become smaller than 
G s tr(A*,P, W), and for sufficiently small k and p near one, 
G(/i, K,p, W) w G str (fj,,p,W), and, as we will show in 
Sec. IV, the simple interpolation formula 

G(n,K,p,W) w max{Goo(/J, K,p),G 3t i(n,p, W)} (2.8) 

provides and excellent description of the simulation data 
with the EMT form for Geo near p = 1 where the finite 
size corrections are the most important. 



B. Kagome lattice energy 

The kagome lattice has three sites per unit cell, which 
we take to be the three sites, labeled 1, 2, and 3, on an 
elemental triangle in Fig. 3. All bonds in the lattice are 
parallel to one of the vectors, ei, e 2 , and e 3 , specifying 
the direction of bonds in the unit elemental triangle. We 
label each site by an index i = (1, a) where 1 = (hjfo) 
labels the position of site I in the unit cell at dimen- 
sionless (i.e., taking a = I) position ri = 2(Ziei + ^2), 
where the factor of two arises because the separation be- 
tween unit cells is twice the bond length. The index 
(7 = 1,2,3 labels the site within the unit cell, and the 
position of site t is = ri + s CT , where si = 0, S2 = ei, 
and s 3 = — e 3 . When the lattice is distorted, maps to 
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FIG. 3. A section of a kagome lattice showing lattice direc- 
tions ei, e2, and e3. Sites and AW bonds in our choice for 
the unit cell labeled 1, 2, and 3 and 1 — 6 in italics, respec- 
tively. Two of the six NNN bonds in a unit cell are shown 
at the bottom of the figure. The other four have a similar 
configurations aligned along e2 and e3 rather than ei. 



a new position = + ui, where is the displace- 
ment vector of site £. The stretching energy of a bond 
(£, £') connecting nearest-neighbor (NN) sites I and £' is 
(l/2)(ju/a)(u«< • i u <) 2 where u u , = u e , - u e and ru> is 
the unit vector parallel to rt> — rg. Bending energy is de- 
termined by angles 6u,i" between NN contiguous bonds 
(£,£') and (£',£") parallel to one of the lattice directions 
e„. To harmonic order, Ou't." = (l/«)(u«' — u^«) ■ e; 
where is the vector perpendicular to e„, which spec- 
ifies the directions of both bonds (£,£') and (£',£"), ac- 
cording the right hand rule about the outward normal to 
the plane. The bending energy arising from a nonzero 
du'l" is (l/2)(fv/a 3 )(#££'£») 2 , and the harmonic energy 
for our latticed is thus 



x 

n ' 



III. EMT: OVERVIEW AND RESULTS 

Effective medium theory [20-22] is a well-established 
approximation for calculating properties such as elec- 
tronic band structure and phonon or magnon dispersions 
of random media. Various formulations of EMT exist, 
but all replace the random medium with a homogeneous 
one whose parameters (such as NN hopping strength or 
spring constant) are determined by a self-consistent equa- 
tion. Here we use the formulation, presented in greater 
detail in App. A, in which self-consistency equations for 
the effective-medium parameters are determined by re- 
quiring that the average multiple-scattering potential or 
T-matrix associated with a single random bond (or group 
of bonds) in the homogeneous effective medium vanishes. 

The development of an EMT for elastic networks with 
both stretching and bending forces presents challenges 
not encountered in networks with stretching forces only. 
In our system, stretching forces are associated with sin- 
gle NN bonds, but bending forces are associated with 
contiguous pairs of NN bonds that couple next-nearest- 
neighbor (NNN) sites in what we call phantom NNN 
bonds that only exist if both its constituent NN bonds 
are occupied. Thus, the replacement of a single NN 
bond, which we will refer to as the replacement bond, in 
an effective medium affects not only the stretching en- 
ergy of that bond but also the bending of the two NNN 
"phantom" bonds containing that bond. There are cur- 
rently two versions [23-26] of EMT on lattices with bend- 
ing energies, which we will refer to as EMT I and EMT 
II, respectively, that deal with this problem in different 
ways. In both methods, the effective medium is charac- 
terized by homogeneous stretching and bending moduli 
fi m and K m respectively. 

In EMT I [24, 25] , the replacement bond has a stretch- 
ing modulus /i s and a bending modulus k s that take on 
respective values /i and K if the bond is occupied and 
if the bond is not occupied. Thus, the probability distri- 
bution for fi s and K s , 

P(fl s ,K s ) = pS(^ s -fl)S(K s -K) + (l-p)S(^J, s )S{K s ) 1 (3.1) 



E — E s + Ef, 

1 9e.e'( u w ' r«') 2 



2 a 
I k 



(2.9a) 
(2.9b) 



BiW = -- V 9t,t'9v,t»{0u'i») 2 (2.9c) 
1,1', I" 

= ^2 91,1' 91', I" [( u «' - uw) ' e^] 2 , 



where gi^ = 1 if the bond connecting sites £ and £' is 
occupied and g^f — otherwise and where the sum in 
Eb(n) is over all continuous triples of sites £,£',£" along 
each filament. The combination fi/a is the stretching 
elastic constant, and the combination n/a 3 is the bending 
constant for contiguous pairs of bonds. 



exhibits strong correlation between the values of /z s and 
k s . As discussed in detail in Ref. [25], the bending con- 
stant of the NNN phantom bonds containing the re- 
placement bond is calculated assuming it is composed of 
two elastic beams connected in series, one with the bend- 
ing modulus K m of the effective medium and one with the 
bending modulus k s of the replacement bond. This leads 
to a bending constant, k c , for both NNN bond contain- 
ing he replacement bond, that is a nonlinear function of 
and Km '. 



1 



1 



(3.2) 



The perturbation to the effective medium arising from 
the replacement bond then consists of a stretching energy 
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on that bond with spring constant (fi s — fj, m )/a and bend- 
ing constants on the two TV TV TV bonds of (k c ~ K m ) /a 3 . It 
turns out, as discussed more fully in Refs. ([24, 25]) and in 
App. B, that the EMT equations in Method I do not close 
unless an additional term, with coupling constant A m , 
coupling the angles on neighboring TV TV TV bonds along a 
single filament is added to the effective medium energy. 
This leads to an additional term in the perturbation aris- 
ing from the replacement bond, with coupling constant 
— X m /a 3 , that couples the two angles on the two NNN 
bonds containing the replacement bond. Thus this EMT 
is characterized by three parameters /i m , re m , and A TO . 

In EMT II [23, 26], the phantom NNN bonds car- 
rying the bending energy are elevated to the status of 
real bonds that exist whether or not the NN bonds of 
which they are composed are occupied. This leads to 
the great simplification that stretching and bending are 
effectively decoupled, and there is no necessity of intro- 
ducing A m . There are separate probability distributions 
for the stretching modulus of the NN replacement bond 
and for the bending modulus of the NNN replacement 
bond. However, since the bending modulus of a NNN 
bond is zero unless both NN bonds comprising it are 
occupied, the probability that the NNN is present with 
a bending modulus k is chosen to be p 2 , the probability 
that any two given TV TV bonds are occupied: 



p6(/j, s - fi) + (1 - p)d(fi a ) (3.3a) 
P 2 5{k s -k) + (1-p 2 )6(k s ). (3.3b) 



The predictions of these two methods are qualitatively 
similar: they both yield a rigidity percolation threshold, 
below which the lattice loses rigidity, and a smoothly 
varying shear modulus approaching the pure lattice value 
of Go as p — » 1. Figure 4 plots values of G(k,n) calcu- 
lated using both EMT methods for different values of 
«;/(/ia 2 ) and p along with the results of simulations. The 
k = simulation curves follows the finite-size result of 
Eq. (2.7). The n/(fJ,a 2 ) = I(T 6 curve follows the EMT 
curve, which corresponds to infinite W, with increasing 
P and then the k = finite size curve after the two cross 
in accord with Eq. (2.8). As is the case for the trian- 
gular lattice [26] , EMT II predicts a value of pb that is 
close to that measured in simulations, whereas EMT I 
predicts a considerably larger value. In addition, for val- 
ues of k/ (na 2 ) beyond 0.1, the EMT I numerical solution 
ceases to exist for small p, and as a result the n/ (/ia 2 ) = 1 
curve is not included for EMT I in Fig. 4. For p > 0.7 
and n/(fia 2 ) > 10~ 6 , the two EMT curves and the sim- 
ulation curves are essentially indistinguishable, and near 
p = 1, they become analytically identical with the effec- 
tive medium modulus \x m satisfying a scaling form 



jU m = /i • *(t) 
where the scaling variable is 



{1-pf 



(3.4) 



(3.5) 




a. 



FIG. 4. Comparison between the shear modulus (normalized 
by the shear modulus Go at p = 1, so it is equivalent to \Xm/ \i) 
obtained from our numerical simulations (data points): EMT 
I (solid lines) and EMT II (dashed lines). Different colors 
mark different values of K,/((ia 2 ), and from top to bottom the 
corresponding values of K,/(fia 2 ) are 1, 10 -2 , 10 -4 , 10 -6 , 0. 
The dotted line indicates the finite size correction predicted 
by Eq. (2.8) As p increases toward 1 the simulation curve 
for n/(^a 2 ) = 10" 6 follows first the EMT curves and then 
the finite-size k — curve beyond the point at which the 
latter curves meet in accord with Eq. (2.8). For the values of 
K/(fj,a 2 ) > 1(F 6 , the EMT solution is larger that the « = 
finite-size result, and thus the simulation data follow the EMT 
curves. 



and the scaling function is 



*(t) = (-1 + VI + At) 2 /{At) . 



(3.6) 



with the constant A = [40(1 - \/2/3)} 2 ^ 53.9. This 
scaling function can be expanded in the following limits 



*(r) 



1 if t > 1 
At/4 if r < 1 



(3.7) 



Therefore, for the case of t 3> 1 corresponding to p near 
1 so that K/(na 2 ) ^> (1 — p) 2 , the shear modulus is ap- 
proximately 



G 



V3 



I'm ^ V3 H 

a 8 a 



Go, 



(3.8) 



indicating that the macroscopic elastic response of the 
network is dominated by the stretching stiffness of the 
filaments and reaches the affine pure lattice limit as p — > 
1, and we shall call this the "stretching dominated" elastic 
regime. 

In the other limit r <C 1 corresponding to smaller k 
or larger 1 — p so that K/(fj,a 2 ) -C (1 — p) 2 , we have the 
shear modulus 



G~ A 



V3 



32 a 3 (l -pf 



A 



V3 k L 2 
32 = 



(3.9) 



indicating that the macroscopic elastic response of the 
network is dominated by the bending stiffness of the 



T 

FIG. 5. (Color online) Scaling behavior near p = 1 for dif- 
ferent values of k = n/(/ia 2 ). The data points for the shear 
modulus stem from our simulations for 200 2 unit cells. The 
continuous curves depict our numerical solutions for /i m of 
the EMT equations for the corresponding values of k. The 
dashed red curve shows the asymptotic scaling function pre- 
dicted by the EMT, see Eq. (3.6). As expected, the simulation 
data and the numerical solution of the EMT equations deviate 
from this scaling function as 1 — p becomes large. 



filaments, and we shall call this the "bending domi- 
nated" elastic regime. The EMT solution for p near one 
is plotted in terms of the scaling variable t together with 
the asymptotic scaling function ^(r) in Fig. 5. Because 
the asymptotic solution (3.6) assumes small (1 — p), it 
requires very small K/(jia 2 ) to make the regime \&(t) ~ 
At I 4 visible. We conclude from this that the elasticity 
of the network is bending dominant as long as 

K/(fia 2 ) < (1 -pf. (3.10) 

As discussed in Sec. I, near the rigidity threshold, 
G(p,k) vanishes as fj,(p — Pb) 1 , where t = 1 and pb — 
0.6920 for EMT I and p b = 0.6180 for EMT II, times a 
scaling function of b m = n/(fia 2 ) that is a constant at 
large b m and proportional to b m at small b m . Thus, at 
sufficiently small b m for all p < 1, response is bending 
dominated with G oc k; at sufficiently large b m , response 
is stretching dominated with G oc \x as shown in the phase 
diagram of Fig. 2. 



IV. NUMERICAL SIMULATIONS 

In the numerical portion of our work, we study the 
elasticity of the filamentous kagome lattice by gener- 
ating diluted lattice conformations on a computer and 
then calculating their mechanical response numerically. 
Practically this is done via deforming the network by 
imposing a certain rj and by then minimizing the elas- 
tic energy (2.9) over the non-affine displacements 5u.g = 
U£ — r\Yg of the sites using a conjugate gradient algo- 
rithm. To explore the response to shear, e.g., we set 



FIG. 6. The bulk modulus B (normalized by the bulk modu- 
lus Bo at p = 1) from our simulations of 48 2 unit cells. 

Vij = l{$ix5jy + SiySj x ), where 7 specifies the magni- 
tude of the applied deformation. We use the same small 
magnitude 7 = 0.01 for all deformations. For a range of 
p-values, we generate up to M = 200 random conforma- 
tions, we calculate several measurable quantities for mul- 
tiple k- values, and we average arithmetically over all con- 
formations. The quantities that we calculate are the elas- 
tic moduli and the corresponding fractions of rigid con- 
formations riijki and non-affinity parameters T^^i ■ Uijki 
is defined as the number of conformations with non-zero 
Gijki (non-zero meaning larger than a small numerical 
threshold, here 10 -8 ) divided by M. The non-affinity 
parameters [40] 

% = ^EW) 2 . ( 4 -i) 

where N is the total number of sites and 8u® is the equi- 
librium non-affine displacement of site I in the presence 
of the r\ that leads to Cijki , measure the deviation from a 
homogeneous strain field. To mitigate boundary effects, 
we apply periodic boundary conditions on all boundaries. 
We simulate system sizes ranging from 6 2 to 200 2 unit 
cells which corresponds to N = 108 and N = 120, 000 
sites, respectively. 

Figures 4 and 6 show log-linear plots of the shear mod- 
ulus G and the bulk modulus B, respectively, as func- 
tions of the occupation probability p for 48 2 unit cells. 
The EMT predicts the bulk modulus to be proportional 
to fj, m , and hence, both /x and B should follow the same 
curves when they are normalized by their respective affine 
values Go and Bo at p — 1. Within the numerical er- 
rors, this prediction is indeed borne out by Figs. 4 (a) 
and (b). Moreover, these figures are consistent with the 
EMT prediction that the rigidity percolation threshold 
Pb be the same for all n > 0. Numerically, we find 
Pb = 0.605 ±0.005 which is only slightly smaller than the 
EMT II prediction pb rs 0.616 and about 15% smaller 
than the EMT I prediction p b rj 0.692. Previous EMT 
predictions for elastic networks have produced somewhat 
larger values for the rigidity threshold than found in cor- 
responding numerical studies, see e.g., Refs. [24, 25], and 
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(a) 




P 

(b) ' 



1 




P 

FIG. 7. (Color online) Non-affinity ratios (a) T xyxy and (b) 
Taxxx for 48 2 unit cells. Color-code is the same as in Figs. 4 
and 6. 



this apparent trend is continued here. For k = 0, the 
EMT predicts a first-order rigidity percolation transition 
at p = 1. The curves for k = in Figs. 4 (a) and (b) rise 
from zero with finite slope below p = 1, a clear finite- 
size effect. Below, we will analyze this finite-size effect 
systematically, and we will find that our data is indeed 
consistent with a first-order transition at p = 1 . 

Figure 7 shows the non- affinity functions T xyxy and 
^xxxx extracted from the simulation runs producing the 
results for \x and B discussed in the previous paragraph. 
For p = 1, the network is affine, and the non-affinity 
functions are zero. The curves for k > are expected to 
have their maxima at pb, and our data is consistent with 
that expectation. The data for n = the curves for G 
and B should show no indication of the existence of pb, 
and indeed they do not. 

Next, we turn to the fractions of rigid conformations. 
These quantities are convenient for the purposes of finite- 
size analysis, in particular for the expected first-order 
transition for k = at p = 1. In the following, we will 
with focus on n = n xyxy . Figure 8(a) displays n for 48 2 
unit cells for n — and several values of k > 0. The 
curves for k > are consistent with the above state- 
ment that the rigidity percolation threshold for all non- 
vanishing values of k is at pb = 0.605±0.005. Figure 8(b) 




0.4 0.5 0.6 0.7 0.8 0.9 1.0 
P 




0.4 0.5 0.6 0.7 0.8 0.9 1.0 
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FIG. 8. (Color online) Fraction of rigid conformations for (a) 
48 2 unit cells and various values of k [with the same color- 
code as in Figs. 4 and 6.], and (b) for k = and system sizes 
of S 2 unit cells with S = 6, 12, 24, 48, 200 (increasing from 
left to right). The bold lines correspond to the theoretical 
result (4.3). 

shows n for k — for a variety of system sizes of S 2 unit 
cells. Along with the data points, the figure follows pre- 
dictions for u(k — 0,p, S) that follow from the elementary 
combinatorics (see Sec. II A). Consider a system such as 
that shown in Fig. 1 with N x = 2S bonds along its two 
sides. Only those filaments, which consist of N x bonds, 
that extend from the top to the bottom of the sample 
will contribute to G, and G will be zero unless at least 
one of the filaments starting on bottom reaches the top, 
which occurs with probability 

P m (p,N x ) = l-(l-p N *) N *. (4.2) 

This function becomes a step function at p=l when 
N x —> co. Thus, we expect the data for n — to be 
fit by the function 

n(K = 0,p,S) = P m (p,2S). (4.3) 

Figure 8 (b) reveals that the data is fit by this predic- 
tion remarkably well even though there is not a single 
adjustable fit-parameter involved. Now, we are in po- 
sition to discuss the infinite-size limit. For S — > oo, 
n(fc = 0,p, S) approaches a unit-step function at p = 1. 
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FIG. 9. (Color online) Scaling behavior of the shear modulus 
near pb- The plot shows both our simulation results and our 
numerical solutions of the EMT I equations relative to the 
asymptotic scaling form near pt, predicted by the EMT, cf. 
Eq. (??). The data points stem from our simulations for 48 2 
unit cells, and we use pb = 0.605, C\ — 0.3, C2 = 0.007, and 
t — 0.2. The lower branches correspond to values of p below 
pt and are a typical finite-size effect. The continuous curves 
represent the numerical solutions to the EMT I equations, 
which reduce to Eqs. (B22) to(B24) as Ap tends to zero with 
ci = 0.03802, c 2 = 4.697, and p b = 0.692. 

This establishes that the rigidity percolation transition 
for k — is a first-order transition at p = 1 where 
the shear modulus jumps discontinuously from zero to 
its affine value. 

Now, let us look at the elastic response near p = 1. 
Figure 5 shows the shear modulus and fj, m as functions 
of the scaling variable r defined in Eq. (3.5). We find 
remarkably good agreement between the simulation data, 
the numerical solution of the EMT equations as well as 
the EMT prediction for the asymptotic scaling function 
M(t) including the predicted value for the constant A. 
Note that this agreement is obtained without adjusting 
any fit-parameters. 

Figure 9, finally, displays the scaling behavior near the 
rigidity percolation threshold pb- Our numerical data col- 
lapses in full qualitative agreement with the EMT scaling 
form (??) albeit with an exponent t = 0.2 that is signifi- 
cantly smaller than the t — 1 predicted by the EMT. The 
value ci = 0.25, on the other hand, is in semi-quantitative 
agreement with the EMT. 

V. DISCUSSION 

We have introduced the diluted kagome lattice in 2- 
dimensions as a model for networks of scmi-flcxible poly- 
mers, whose maximum coordination number is four. We 
identify straight lines of contiguous occupied bonds as 
filaments that, following previous work [7—11, 19], we en- 
dow with the energy of an elastic beam characterized by a 
stretching (Young's) modulus fi and a bending modulus 



k. We contrast this model, which has filaments of ar- 
bitrary length, with the Mikado lattice whose filaments 
are of finite (and usually fixed length). We show that 
when k = 0, this kagome model has affine response and 
non- vanishing elastic moduli so long as there are sample- 
traversing filaments, which is the case for samples with 
finite lengths and widths W, even for bond-occupation 
probabilities p less than one. As W increases, the prob- 
ability of sample-traversing filaments decreases, and in 
the W — > oo limit, elastic moduli, which are nonzero 
for undiluted p = 1 lattice, fall precipitously to zero for 
p = 1~ in a first-order transition. The undiluted lattice 
has ^-independent, and thus affine, elastic moduli. The 
addition of bending forces restores rigidity for any n for 
p greater than pb, the rigidity percolation threshold, and 
elastic moduli approach the non-zero affine values of the 
p = 1 lattice as p —> 1. We argue that this is the un- 
derlying cause for the affine limit found in the large l/l c 
limit (p near one) in the Mikado model. 

We use two recently introduced lattice-based effective 
medium theories to calculate elastic moduli of the diluted 
kagome lattice as a function of fx, k, and p, we we calcu- 
late scaling forms for the shear modulus both near p = 1 
andp = pb- Both forms are a function of the unitless ratio 
b m = K / (M a2 ): where a is the lattice spacing, and yield a 
crossover for all pb < p < 1 between bending dominated 
response at small b m and stretching dominated response 
at large b m . We supplement our EMTs with numerical 
simulations of the shear modulus and other functions, 
such as that measuring the degree of nonaffine response. 
The results of these simulations agree qualitatively with 
EMT predictions for all p and quantitatively with them 
near p = 1. 

Our study of the kagome lattice provides insight into 
the behavior of 3-dimensional filamentous lattices with 
maximum coordination number of four. With stretching 
forces only, these lattices are subisostatic with of order 
one zero mode per site, and one might, therefore, conjec- 
ture that their elastic moduli should vanish when k = 0. 
However, simulations [33, 34] on two model lattices con- 
sisting of straight filaments yields curves of elastic moudli 
as a function of p that look very similar to those of the 
kagome lattice presented here with an approach to affine 
k- independent response as p — » 1. The underlying cause 
of this result is that straight sample-traversing filaments 
in any dimension give rise to nonvanishing macroscopic 
elastic moduli. Indeed, exact calculations [33] for one of 
the 3<i-lattices with elastic-beam energies show that all 
elastic moduli are nonzero when n = and p = 1 and 
that response is affine. Nearp = 1, simulations on both of 
the 3d lattice show a crossover from a bending dominated 
regime with G ~ (K/a 3 )(L/a) 2 to a stretching dominated 
regime in which G ~ /i in agreement with the predictions 
of our kagome-based lattice EMT. In fact the results of 
one of the simulations [33] track the scaling function of 
the our kagome EMT. We conjecture that the impor- 
tant common feature of the 2d and 3c? lattices is that 
they both consist of sample-traversing straight filaments 
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in the undiluted limit. Away from p = 1, bending forces 
stabilize the lattices in a process that is not so sensitive to 
lattice dimension, even though the central-force 3d lattice 
is subisostatic. Our EMTs provide some support for this 
conjecture. The characteristic property of the kagome 
EMT solutions, such as the form of the scaling function 
near p = 1, are a direct consequence of the existence 
of lines of zeros in the phonon spectrum at k — that 
get raised to nonzero frequency for k > as shown in 
Fig. 10. This feature causes the integral 7?2,s [Eqs. (B14) 
to (B17)] to diverge as \fb~^ x as b m — n/(pia 2 ) — > 0. In 
three dimensions, the lines of zero modes in the phonon 
spectrum become planes of zero-modes, and we believe 
that the 3d version of i?2,s will have a from analogous to 
Eq. (B17) but with g(q y ) replaced by a similar function 
of q y and q z : 



H 2 ,s{bm,0)~ / d 3 q 



\{B>f)l 



cql + b m g{q y ,q z ) 



-1/2 



(5.1) 



When b m — 0, the denominator vanishes for all q y and 
q z when q c = 0, and the integral diverges. Nonzero b m 

— 1/2 

yields a b m divergence just as in the 2d kagome case. 
Unfortunately, the 3d lattices are quite complicated, with 
a 54-site unit cell in one case [33] and complicated cross- 
ing configurations in the other [34] , and we have not been 
able to set up a 3d EMT to verify this conjecture. 

Much of the work on the elastic properties of filamen- 
tous networks (with the notable exceptions of [9, 32, 41]) 
have focused on models, such as the Mikado or kagome 
lattice presented her, consisting of straight filaments. 
This work has also largely focused on what are really 
mechanical models with beam energies assigned to the 
filaments and the effects of thermal fluctuations ignored. 
As our analysis shows, these are exceptional lattices be- 
cause they can support stress with central forces only 
even when their coordination number z is substantially 
less the Maxwell stability limit of z = 2d. Filaments in 
real biological gels are not straight, and their energies 
are not described by that of an elastic beam. It is thus 
legitimate to ask how seriously one should take models 
such as those presented here as realistic descriptions of 
real bio-gels of semiflexible polymers. If constituent fila- 
ments are not straight, then in general, at least some of 
elastic moduli of even perfect lattices vanish in the ab- 
sence of bending forces. Prime examples of this behavior 
are found in the two-dimensional honeycomb lattice and 
the three-dimensional diamond lattice whose bulk moduli 
are nonzero but whose shear moduli vanish in the n — > 
limit, but it has also been observed in a more realistic 
model of a filamentous lattice [41]. If k is needed to sta- 
bilize the system, the elastic response to external stresses 
will involve bending and thus be nonaffine. So, it would 
seem that the straight-filament models are not such good 
ones for real systems, though they do provide us with 
valuable insight into how network architecture influences 
clastic response. One of the predictions of the straight- 



filament models is that response becomes more stretch- 
ing dominated and more affine as k increases. Thus, it 
is plausible that if k is large enough, the elastic response 
of even those lattices whose shear moduli vanish if k = 
will exhibit stretching dominated, nearly affine response 
for sufficiently large n. 
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Appendix A: EMT Generalities and the EM 
Dynamical Matrix 

To implement the T-matrix version of EMT, we deal 
with the dynamical matrix D m of the EM, its associated 
Green's funcition G m = — (D m ) _1 (we consider only zero 
frequency), the perturbation V associated with bond re- 
placement, the dynamical matrix D = D m + V, and 
its associated Green's function G = — D _1 . [For more 
details of this formalism, see Refs. [42] and [25]] Because 
there are three sites per unit cell, all of these matrices are 
6N x 6N matrices, to be detailed further below, where N 
is the number of unit cells in the lattice. We will specify 
G m more completely below. With these definitions, 



G = G m + G m • T • G m = (G m - V)" 



where 



T = V- (I-G m • V)" 



(Al) 



(A2) 



The EMT self-consistency equation requires that the 
disorder average over the probability distribution of 
Eq. (3.1) for EMT I or Eq. (3.3a) or EMT II of the per- 
turbation vanishes 



(T(/i„« a )) =0. 



(A3) 



All of the matrices discussed here can be expressed in 
a position or wavenumber representation. The six inde- 
pendent displacements in a unit cell can be expressed as 
a six-dimensional vector Ui = (111,1,111,2,111,3) for each 
1 or as its Fourier transform U q = ^ e~ lq r<! U(l) 



(u q ,i,u qi2 ,u qi3 ) 
tem is then 



The energy of the bond-replaced sys- 



£? = ^U,D Iir U{ 



i.i' 



i^U_ q D q , q ,U q ,, 



2N 2 



(A4) 



11 



where D\_y and D qq / are 6x6 dimensional matrices for 
each pair (1, 1') and (q, q'). The effective medium is trans- 
lationally invariant and 

D™ q ,=M q , q ,D™. (A5) 



The energy of the effective medium is constructed by 
occupying all bonds with identical beams with stretching 
and bending moduli [i m and n m and adding (for EMT I) 
an additional coupling of strength X m coupling angles in 
neighboring NNN bonds along a filaments [See Ref. [25] 
for details]. The EM energy is then 



E m = E s (^ n ) + E h {n m ) + E%(\ m ) (A6) 

^E u q - D r u - q ( A7 ) 

q 



where E s (/j, m ) and Eb(n m ) are evaluated at gi t i> = 1 for 
all bonds (£,£'), where 



E bb( X m) - -f" ^lA,<3^,<3,<4' ( A8 ) 

where it is understood that £i , f 2 , £3 and £4 are all con- 
tiguous sites along a single filament. G™ can be con- 
structed from the Fourier transforms of the stretching 
and bending energies on NN and phantom NNN bonds. 
For example, the stretching energy of bond 5 in Fig. 3, 
connecting site £ = (1, 3) at position ri — e 3 with site 
£' = (1, 2) at position ri + 2e2 + ei such that r' e — ri = e 2 
is |(/U m /a)[(e^ — e^) • e 2 ] 2 . The Fourier transform of 
(e £ / - ei) ■ e 2 is 

e - iq .l (e - iq .e 2u2 q _ U3 q) . e2 = e -iq r 1B 2 ^ . ^ (AQ) 

where B§ q is specified in detail in Eq. (All) below. Simi- 
lar procedures apply to all stretching and bending bonds, 
and the EM dynamical matrix can be expressed as 

6 

- — 2^ B ™q B n,-q 

n=l 
6 

I K ™ \^ R 6 R b 
a 3 A^i ■ D n,q J -'n,-q 

n=l 

+ ^E 2C0S (q- e «) B «,q B «.-q' ( A1Q ) 



where 

B 5 liq = {-ei,ei,0,0} 

B|, q = {0,0,-e 2 ,e 2 } 

B^ q = {e 3 ,0,0,-e 3 } 

BI iq = {e- iqei e 1 ,-e 1 ,0,0} 

B;? jq = {0,0, e -^ e2 e 2 ,-e 2 } 

B 6, q = {-e 3 ,0,0,e-^ e3 e 3 } (All) 

and 

B^ q = 2{2e r ,-(l + e ^ ei )e r ,0,0}, 

B^ q = 2{0,0,2e 2 L ,-(l + e it »- e2 )e 2 L }, 

B^ q = 2{-(l + e ^ e3 )e 3 L ,0,0,2e 3 L )}, 

B 4, q = 2{-(l + e- jqei )e 1 ,-2e 1 ,0,0}, 

B| iq = 2{0,0,-(l + e - it »- e2 )e 2 L ,2e 2 L )}, 

B^ q - 2{2e^ 0, 0, -(1 + e~^)ei)}, (A12) 

where are the unit vectors perpendicular to the bonds, 
e„ should be understood as {e„ i2: , e n , y }, and the same for 
so all these vectors are 6 dimensional. The factor of 
2 is from the fact that the length of the bonds is 1/2. 

The scattering V can also be expressed in this form. 
We assume that the changed bond is bond 1 in the unit 
cell at r = as marked in Fig. 3. Then we have 

\T _ Ms ~ Mm ps r, s 

+ ~3 l B l,q B l,-q + B 4,qB 4 ,- q J 

-^(B^X-q + BtqB?,-,) (A13) 

Appendix B: Asymptotic solutions of the EMT I 
equations 

Asymptotic solutions for the EMT I self-consistency 
equation (A3) has been developed in Ref. [25] in the limit 
of small k/ jj,. In this section we review this asymptotic 
solution and apply it to the case of kagome lattice. For 
convenience we define the notation 

bm = K m/ (Pm®?) 
^m — ^m/ (Mm^ ) 

b = K/(/j, m a 2 ) 

(Bl) 

and 

H (b m , l m ) — ~ G (/imi ^m) 

= -G(l,b m ,l m ), (B2) 

where the second line is derived from the definition of the 
dynamical matrix. The self-consistency equation (A3) 
can be solved by projection to the space |Bf), \B\), |B|) 
that spans V. In this basis we can rewrite the EMT 
matrix equation into three independent equations 
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Mm = M" 



p- Hi (b m , l r , 



1 — Hi {b m , l m ) 

- l m - 2( ^ + -!-) 1 (Z ro tf 2 + &m# 3 ) + 2b m l m H 2 + (b 2 m + l 2 m )H 3 = 0, 



(B3a) 
(B3b) 
(B3c) 



where i?i and #3 are the projections of the Green's func- 
tion defined as 



ffi(M) 



(b;|h(&,o|b;> 
(b?|h(&,o|b5) 
(b5|h(6,o|b5), 



(B4) 



with the inner product defined as 
<B?|H( M )|B?) 

=^E B U- H ( & '0-Bt c 

q 

=TrH(6,Z)[|B?)(B?|] 



(B5) 



where the sum is over all N vectors in the first Brillouin 
zone of the kagome lattice, and the trace is understood to 
include the sum over these vectors (along with the factor 
of l/N in addition to the sum over the 6 dimensional 
space of Bf . These equations are exactly equivalent 
to the matrix equation (A3). 

For the special case of k = 0, equations (B3a, B3b, 
B3c) simplify and give 



Mm = M 

b m = 
lm =0 



, p-gi(0,0) 
1 - Hi (0, 0) 



(B6) 

(B7) 
(B8) 



As discussed in Ref. [25], it is straightforward from the 
definition D and the fact that the central force undiluted 
kagome lattice is isostatic with z — 2d to derive the rela- 
tion 



ffi(0,0) = l, 



(B9) 



The effective medium filament stretching stiffness is then 
Mm = for p < 1 and jj, m = [i for p = 1. Therefore 
this k — EMT solution indicates a first order rigidity 
transition at p = 1. 

In the following we solve these equations asymptoti- 
cally at small k near the two critical points p = 1 and 
Pb- 



1. Asymptotic solution near p = 1 

The EMT solution for small k > can be calculated 
perturbatively from the k = solution (B6). Near p = 1, 



as discussed in Ref. [25], we can make simplifications to 
the self-consistency equation (B3b,B3c) using the fact 
that k< 1 so that 



K m = n(2p- 1) 

Am = 0, 



(B10) 



(which we shall verify later) and therefore we only need 
to solve Eq. (B3a) using perturbation. 

In contrast to the perturbative calculation in the tri- 
angular lattice, the kagome lattice effective medium is 
isostatic as n m -4- 0, and thus the phonon Green's func- 
tions exhibit singularities. These singularities correspond 
to the zero- frequency floppy modes of the dynamical ma- 
trix, and make diverging contributions to H. 

Nevertheless, perturbation theory around the k = 
solution is still well defined. We shall see below that the 

— 1/2 

H singularity is proportional to b m and thus all the 
terms in the self-consistency equations are non-singular. 
The expansion of H\ at small b m can be calculated using 
the following equality 

<B;|H(6 ro , 0)|BJ) + 6 m (B?|H(6 m , 0)|Bj) 
= TrH(6 m ,0)(|B^)(BJ| + 6 m |B^)(B?|) 



= -TrH(6 m ,0) 



^(|B^)(B^|+6 m |B^)(B^|) 



.71=1 



(Bll) 



These relations follow because all of the 6 NN and the 
6 NNN bonds in a cell are, respectively, equivalent by 
symmetry allowing the trace over one set of bonds to 
be replaced by 1/6 the trace over the sum of the bonds. 
But the quantity in square brackets is just the inverse of 
H(6 m ,0), and the final result follows. 

Employing the analysis of the phonon modes in 
Ref. [42], we find that |B*) has a nonzero projection onto 
the floppy mode branch, which we call \ vf), whereas |Bf) 
docs not. The floppy modes branch has low frequencies 
that are proportional to \fb~m~ along symmetry directions 
at {f , f , ^, ^, ^F, in the first Brillouin zone, as 
shown in Fig. 10. In this calculation we shall use the 
direction n/2 which correspond to q x = as an example. 
Near the q x = line the floppy modes branch phonon 
Green's function takes the form, 



Gfl. 



1 



oppy,§ ,q 



T&Ql + 16b m g(q y ) ' 



(B12) 
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where 



9(ly) 



3Qm — 2Qs — 1y 

Qm — Qs 



(B13) 



with Qm = 2tt/s/3 and Qs — 4(36 m /2) 1 / 4 and thus g(q y ) 
takes value between 2 and 3 in the first Brillouin zone. 
This form of the Green's function is shown in Fig. 10. 
It was derived in Ref. [42] in which the weak additional 
interactions are NNN bonds rather than bending forces. 

It is straightforward to calculate the singular part 
of (B*|H(b m ,0)|B5) along symmetry directions of the 
floppy modes, which involves the integral 



H 2 ± a (b m ,0) 



= *'() 



d 2 q |(Bj>/)| 2 
lBZ,\q x \<\q y \/V5 ( 2?r ) 2 Telx + b mg(q y ) 



(B14) 



where Vq — V3/2 is the area of the kagome lattice unit 
cell, and the condition \q x \ < \q y \/\/3 confines the inte- 
gral to be around the floppy mode directions ir/2 and 
37r/2 near which the expansion (B12) is valid. The total 
contribution involves the singular part of all the isostatic 
directions 



H 



2,8 



H 2 



Ho 



Ho 



(B15) 



All these terms exhibit similar behavior so we just discuss 
i?2,f ,s for convenience. In Eq. (B14), assuming b m <C 
1, one can first integrate out q x using contour integrals, 
which gives 



Ho 



1 



(2^) 2 /«o Jo 

6 -l/2 



% dq !ZL !KaMI! 

V V3 y / b m g(q y ) 

(B16) 



— 1/2 

indicating a leading order divergence b m of H 2 at small 
K m . A calculation including all terms in Eq. (B15) yields 



H 2 , s (b m ,0) 



VA b -i/2 
2 m ' 



where 



From this we get 

ili(6 ro ,0)~l 



/ 2/3)- 



51/2 



(B17) 



(B18) 



(B19) 



at small n m . Plugging this back into Eq. (B3a) , we arrive 
at the solution 



H A(2p- l)K/(fj,a 2 ) 




A(2p~ l)n/(p,a 2 

(ApJ 2 



(B20) 



where Ap = 1 — Ap and this solution is asymptotically 
accurate for n <C 1 and Ap <C 1. 




(a) 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 

q y 

(b) 




FIG. 10. (a) Dispersion relation of the floppy mode taking 
jU m = 1, Km = 0.0005, A m = 0, where the frequency uja is 
the square root of the lowest eigenvalue of D(/i m , K m A m ). (b) 
The floppy mode branch of the eigenvalues of the dynamical 
matrix along the direction of n/2 (which is q x = 0) taking 
|U m = l,«m = 10 _5 ,A m = 0. The red solid curve shows 
the actual eigenvalues calculated from the dynamical matrix 
numerically, and the blue dashed line shows the approxima- 
tion we used in Eq. (B12) that uj\ — qy/16 if q y < Qs and 
u)\ — 16b m g{q y ) otherwise (only the q y > Q s part is included 
in Eq. B12 because other contributions are much smaller), 
(c) Comparison between numerical (blue points connected by 
line) and the asymptotic form (B19) (red line) of 1 — Hi(b m , 0). 



2. Asymptotic solution near p = pb 

The rigidity threshold pb can be obtained by plugging 
fi m = into Eqs. (B3a,B3b,B3c), which leads to 

p b ^ 0.6920, 
b m ~ 0.03788, 

l m ~ 0.008032. (B21) 
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Following the calculation in Ref. [25] we arrive at the 
asymptotic solution for fi m 



fj, m = iu,$(K/(pa 2 ))(p - p b ) 



where 



with 



c 2 x 
c\ + x 



(B22) 
(B23) 



ci ~ 0.03802 
c 2 ~ 4.697. 



(B24) 



Appendix C: Asymptotic solutions to EMT II 
equations 

The self-consistency equation of the EMT II can be 
written as [23, 26] 



M 

Km 
K 



p — a* 
1 - a*' 
p 2 -b* 
1 - b* 



(CI) 



where the variables a* and b* correspond to the integrals 
we defined as 



a* =Hi 

b* = (K m /^ m )H 2 



(C2) 



as shown in Ref. [25]. 

Using the asymptotic forms we obtained in 
Eqs. (B17,B19), we have 



p-l + iVAWiKm/flm) 1 /' 



K m p 
K 



- (VA/2)( Km /^„Q 1 / 2 
1 - (V^/2)(« m / Mm ) 1 /2 



(C3) 
(C4) 



Close to p = 1, if we take n/p, <C 1, then 
(VA/2)(n m / Mm) 1/2 < 1 and we have 



Km/ K — 1- 



(C5) 



Plugging this back into Eq. (C3) we get 



Mm = P - 1 + (VA/2)(k/ (l m ) 
M ' ~ (V^/2)(«/ Mm )i/2 

1/2 



1/2 



(C6) 



which is a quadratic equation in [im , and the solu- 
tion to this leads to the same asymptotic solution as in 
Eqs. (B20.3.4). 

Next we discuss asymptotic behaviors near the rigidity 
threshold pi, t 2 (the subscript 2 is used to distinguish it 
from the threshold pb m EMT I). The value of pi, 2 in 
EMT II can be readily obtained from 



0=p-a*, 
0=P 2 -b*, 



(C7) 



along with the condition that a* +b* = 2d/z, as discussed 
in the App. D of Ref. [25]. We then have 



Pb,2 



0.6180. 



(C8) 



The EMT II self-consistency equations CI can then be 
expanded at small 



and this leads to 



with 



Ap =p-Pb,2, 



H<S>2 (n/i/ia 2 )) {p-pba) 



<5>2(x) 



(2p b ,2 + l)x 
Pb,2 b b,2 + (1 ~Pb,2)x 



(C9) 



(C10) 



(Cll) 



This is in the same scaling form as in EMT I [Eqs. (??) 
and (B23)], apart from different constant factors. 
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